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Abstract 



Within this paper we outline a method able to generate truly minimal 
basis sets which describe either a group of bands, a band, or even 
just the occupied part of a band accurately. These basis sets are the 
so-called iVMTOs, Muffin Tin Orbitals of order N. For an isolated 
set of bands, symmetrical orthonormalization of the iVMTOs yields 
a set of Wannier functions which are atom-centered and localized by 
construction. They are not necessarily maximally localized, but may 
be transformed into those Wannier functions. For bands which over- 
lap others, Wannier-like functions can be generated. It is shown that 
TVMTOs give a chemical understanding of an extended system. In 
particular, orbitals for the it and a bands in an insulator, boron ni- 
tride, and a semi-metal, graphite, will be considered. In addition, we 
illustrate that it is possible to obtain Wannier-like functions for only 
the occupied states in a metallic system by generating iVMTOs for 
cesium. Finally, we visualize the pressure-induced s — > d transition. 

Keywords: Wannier functions, density functional calculations, elec- 
tronic structure, muffin-tin orbitals, solid-state 
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1 Introduction 



The electronic structure of condensed matter is usually described in terms 
of one-electron basis sets. Basis functions used for computation are often 
simple, e.g. Gaussians or plane waves, but their number is large, 1-2 or- 
ders of magnitude larger than the number of electrons to be described. This 
is so because the potential in the effective one-electron Schrodinger equa- 
tion, arising say in density-functional theory, is rather deep inside the atoms. 
The results of this kind of computation therefore require interpretation and 
simplification in terms of a small set of intelligible orbitals. The results of 
band-structure calculations for crystals, for instance, are often parametrized 
in terms of tight-binding models. 

The Car-Parinello technique for performing ab initio molecular-dynamics 
simulations using density-functional theory, pseudopotentials, plane-wave ba- 
sis sets, and supercells, [1] created a need to visualize chemical bonds, and 
this caused renewed interest in generating localized Wannier functions for the 
occupied bands. For a set of M isolated energy bands, Sj (k) , the Wannier 
functions, w m (r — T) , enumerated by m (= 1, .., M) and by the N (— > oo) 
lattice translations, T, is a set of orthonormal functions which spans the 
eigenfunctions, 

M 

* jk (r) = J\T*££ w rn (r - T) u jm exp ik • T (1) 

T rn=l 

with eigenvalues Sj (k) , of the one-electron Schrodinger equation. The Bloch 
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states are delocalized, and the set is enumerated by the band index j and 
the Bloch vector k. Wannier functions are not unique, because performing 
a unitary transformation, W m m i. T _ T >, of one set produces another set which 
also satisfies Eq. (fl)l. merely with different jk-dependent phases of the Bloch 
functions. For molecules, it had long been recognized that chemical bonds 
should be associated with those linear combinations of the occupied molecular 
orbitals which are most localized, because those linear combinations are most 
invariant to the surroundings. [2] For infinite periodic systems, Mazari and 
Vanderbilt have developed a useful method for projecting from the Bloch 
states a set of so-called maximally localized Wannier functions. [3-5] 

This paper deals with a different kind of basis set, specifically, mini- 
mal basis sets of localized orbitals, the newly developed iVMTOs (Muffin 
Tin Orbitals of order N, also known as 3 rd generation MTOs. [6-9]) We 
shall demonstrate that with iVMTOs it is possible to generate Wannier func- 
tions directly, instead of via projection from the delocalized Bloch states. 
iVMTOs are constructed from the partial-wave solutions of Schrodingers 
equation for spherical potential wells (overlapping muffin tins) and iVMTO 
sets are therefore selective in energy. As a consequence, one can construct 
an iVMTO set which picks a specific set of isolated energy bands. Since 
iVMTOs are atom-centered and localized by construction, they do -after 
symmetrical orthonormalization- form a set of localized Wannier functions 
which, if needed, can be recombined locally to have maximal localization. 
The iVMTO technique is primarily for generating a localized, minimal ba- 
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sis set with specific orbital characters, and it can therefore be used also to 
pick a set of bands which overlap other bands outside the energy region of 
interest. The corresponding iVMTOs -orthonormalized or not- we refer to 
as Wannier-like. 

iVMTO-generated Wannier functions have so-far been used only in a few 
cases to visualize chemical bonding [9,10] and more often to construct Hub- 
bard Hamiltonians for strongly correlated 3<i-electron systems. [10-13] In the 
future, it may be possible to use iVMTOs for real-space electronic-structure 
methods in which the computational efford increases merely linearly with 
system size, so-called order- N methods. [14, 15] 

In this paper we will focus mainly on showing how iVMTOs may be 
used to pick specific states in insulating, semi-metallic and metallic systems, 
namely boron, graphite, and cesium. Since iVMTOs are unfamiliar to most 
readers of this issue, and more complicated than plane waves, we start out by 
illustrating the main ideas by performing an elementary analytical calculation 
of the 7r-bonds in the simplest tight-binding (TB) model of benzene. Then 
follows a concise summary of the iVMTO formalism. 

2 ATMTO basics 

2.1 Tight-binding calculation of the benzene 7r-bond 

The simplest TB model for the six 7r-electrons in benzene has six orthonormal 
p z -orbitals, <fi, ip^, placed on the consecutive corners of the hexagon. The 
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hopping integrals over the short (double) bonds are: H 12 = H u = if 56 = 

— (1 + d) , and those over the long (single) bonds are: Hqi = H23 = = 

— (1 — d), where d is the dimerization. 

In order to calculate the Wannier function for the three occupied bonding 
states, it is convenient to partition the orbitals into those on the even (e)- 
and those on the odd (o)-numbered sites. The eigenvalue equations are then: 
H 00 u + H oe u e = el 00 u and H eo u Q + H ee u e = el ee u e , in terms of the 3x3 
blocks of the Hamiltonian and the unit matrices. Solving the last set of 
equations for the eigenvector for the even orbitals yields: 



and inserting in the first equations results in the eigenvalue equation: 



for the (Lowdin) downfolded Hamiltonian. 

In the present case where there is no hopping between even or odd sites, 
H ee = H OQ = 0. Moreover, 



(e — H ee ) 1 H eo u, 



(2) 



[H 00 + Hoe (^lee ~ H ee ) 1 H e0 \ U Q = El 00 U, 



(3) 



(1 + d) 







\ 



eo 








(1 + d) -(1-d) 
-(1 + d) j 
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and 



H 00 + H oe (elee — H ee ) 1 if e 



2 \ 



2 (i + rf 2 ) i - d 2 i - d 

1-d 2 2(l+d 2 ) 1-d 2 

y i-d 2 i-d 2 2(i + d 2 ) ) 



(4) 



The latter, downfolded Hamiltonian is periodic with period 3 and is therefore 
diagonalized by the unitary transformation 



U, 



ok 



( X. =1 =±\ 

V3 V2 V6 

11-1 

V3 V2 Ve 

i n 2 

V 73 76 / 



(5) 



yielding singly degenerate a-states (k = 0) with e = ±2 and doubly degen- 
erate e-states (k = ±1) with £ = ±y/l + 3c? 2 . The even components of the 
eigenvectors are obtained from equation ((21): U e f. = e~ 1 H eo U k, and finally 
we can renormalize: u = U / a/2. 

Having found all six Bloch eigenstates, we need to form three Wannier 
functions, that is, three congruent linear combinations of the bonding states. 
The three downfolded p z -orbitals, Xo, defined by: XoU k = foU k + feU e k, 
are obviously congruent. Moreover, they are localized in the sense that xi 
vanishes on the other odd sites (3 and 5). In the present case, they are 
also orthonormal and, hence, Wannier functions. Left-multiplication with 
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(U ok ) 1 = U ko yields: 

r /9 \ /9 \ i 1 

(6) 



1 

Xi 



V2 



V 9 ! + ( \ + d J ^2 + ( | - d ) <^6 - 77^4 
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and X3 and Xs by cyclic permutation of site indices. Here and in the following, 
we work merely to first order in the dimerization, d. The Wannier function in 
Eq. (J1)J) is essentially the iVMTO. It is atom-centered and, as the dimerization 
increases, it becomes lopsided towards site 2, i.e., it spills into the short bond. 
It breaks the symmetry (when d ^ 0) because it was chosen to vanish on the 
odd sites different from its own, and it is not maximally localized (unless d = 
0). Nevertheless, it is a fairly simple matter to achieve maximal localization 
and, hence, to restore the symmetry, by finding a unitary transformation, 
Wt-t 1 , which maximizes e.g. (\w\ 4 ) = Yfn=i \ w t ■ Here, w\ = ^2 t Wi-tXt 
is the maximally localized Wannier function. In the present case, W has only 
one independent matrix element and we find: 



L/3 + 1 + ^-dj ( Vl + ip 2 ) + (l - V^rf) (<^3 + 



~(V3- l-^^ 4 + ^ 5 ) J, (7) 

which is clearly symmetric (bond-centered). From (jHJ): (|x| 4 ) — if + O (d 2 ) , 
and from (jTJ): (|w| 4 ) = l^+^^-d+O (d 3 ) . Hence, for d — 0, the atom-centered 
and the bond-centered Wannier functions are both maximally localized and 
symmetric. 
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Now, the iVMTO set, x° > is obtained without solving the eigenvalue 
equations, i.e., is not obtained by projection from the Bloch states through 
multiplication by (t/fc ) _1 , but in a more tricky way: We first define a set of 
downfolded, energy-dependent orbitals, 

(fio (e) = ip l 00 + tp e (sl ee ~ H ee )~ l H eo , (8) 

which are localized when e does not coincide with an eigenvalue of H ee . 

4> (e) ) = oo , so we realize 



Projection onto the even sites yields: yp e H — e 
that the functions of the set <ft (e) are solutions of the impurity problems 
specified by the boundary conditions that (f> (e) vanishes at the other odd 
sites and is normalized to ip at its own site. Projection onto the odd sites 
yields: 



H>0 



H-e 



H 00 +H oe (el ee — H ee ) 1 H eo — el OQ = — G oa (e) , (9) 



and comparison with (jSJ) shows that, if e is an eigenvalue of the (downfolded) 
Hamiltonian and u Q an eigenvector, then O (e) u Q is an eigenf unction. G (e) 
defined in Q is the resolvent. Finally, we need to find an energy- independent, 
iVth-order approximation to the set 4> (e): We form the set, (fi (e) G OQ (e) , 
of (contracted Greens) functions and add an analytical function of energy 
determined in such a way that the two sets of functions, <p (e) G oa (e) and 
Xo N ^G 00 (e) , coincide when e is on an energy mesh, e , ■••ejv, specifying the 
energy range of interest. By taking the highest-order finite difference on this 
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mesh, one obtains: 



A o 



A N G 00 A N 4> o G 



A[Q...N] A[0...N]' 
which determines the set of (non-orthonormal) iVMTOs. For N 
instance: 



(10) 
1, for 



= [<f> (e x ) G 00 ( ei ) - O (e ) G 00 (e )] [G 00 (e 1 ) - G 00 (eo)p 1 . (11) 



For the simple benzene model, Eq. (jHj) yields: 



, , v 1+d 1-d 
J i (e) = Vi ^2 y?e, 



and 03 (e) and 5 (e) by cyclic permutation of site indices. To order d, the 
downfolded Hamiltonian (Jl| is independent of d, and by subtracting e and 
inverting, we find: 



( e 2 - 3 1 



Goo (£■) 



(e 2 - 1) (e 2 - 4) 



1 \ 

1 e 2 - 3 1 
1 1 e 2 - 3 



Specializing to iV = 1, the LMTO found from (JTTJ) is: 



(i) 

Xi = <Pi 



eo + ei 



( £o e 1+ 4)( £o£l +2 + (£o£l+ 4)^ 2 + 

[e £i + 2 - (e ei + 4) d] ^ 6 - 2^ 4 }, 
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and if, with the benefit of hindsight, we choose e = —2 and e\ = — 1, 
we obtain the exact result (JUJ), apart from the normalization, l/y/2. With 
other choices, the LMTO set is an approximation to the exact Hilbert space 
spanned by Xo or w , as explained in connection with Eq. (fT3j) below. We 
leave it to the reader to convince himself that had we chosen €q = 2 and 
e\ — 1, we would have obtained the Wannier functions for the antibonding 
levels. 

By considering a simple TB model, we have thus learned that the 7VMTO 
procedure for constructing a minimal basis set, specifically a set of localized 
Wannier functions, consists of the following steps: 1) Downfolding to a small 
set of energy- dependent orbitals and 2) a polynomial approximation of the 
latter. The resulting iVMTO set is not orthonormal in general, but may be 
symmetrically (Lowdin) orthonormalized in a third step. Wannier functions 
which are maximally localized, and therefore not symmetry breaking, may 
be obtained in a fourth step. None of these steps require knowledge of the ex- 
tended (Bloch) eigenstates. Although of utmost importance for applications, 
steps (3) and (4) are not specific for the ./VMTO method and fairly standard. 
As a consequence, they will not be considered further in this paper. 

2.2 iVMTOs for real systems 

For real systems, the iVMTO method constructs a set of atom-centered 
local-orbital basis functions which span the solutions of the one-electron 
Schrodinger equation for a local potential, written as a superposition, 
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^2 r vr (|r — R|), of spherically symmetric, short-ranged potential wells, a 
so-called overlapping muffin-tin potential. This is done by first solving the 
radial Schrodinger (or Dirac) equations numerically to find ipm (e, |r — R|) 
for all angular momenta, I, with non-vanishing phase-shifts, for all potential 
wells, R, and for the chosen set of energies, e = eo, cjv. 

The partial-wave channels, Rim, are partitioned into active (odd, in the 
benzene example) and passive (even). The active ones are those for which 
one wants to have orbitals in the basis set. E.g. for the red 7r-bands in Figure 
n the active channels are p z on all the carbon atoms, whereas for the black 
bands, the active channels include all nine s, p, and <i-channels on all the 
carbon atoms. 

For each active channel, Rim, a kinked partial wave, (pjn^ (e, r) (Eq. (jHJ) 
in TB theory), is now constructed from all the partial waves, 



ifRi (e, |r — R|) Y\ m \ y — RJ , inside the potential-spheres, and from one solu- 
tion, ipRjm {e, r) , of the wave-equation in the interstitial, a so-called screened 
spherical wave. The construction is such that the kinked partial wave is a so- 
lution of Schrodinger's equation at energy e in all space, except at some hard 
screening-spheres -which are concentric with the potential-spheres, but have 
no overlap- where it is allowed to have radial kinks in the active channels. 
In passive channels, the matching is smooth. 

It is now clear that if we can form a linear combination of such kinked 
partial waves with the property that all kinks cancel, we have found a solu- 
tion of Schrodinger's equation with energy e. In fact, this kink-cancellation 
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condition (Eq.s (JHJ) and © in TB theory) leads to the classical method of 
Korringa, Kohn, and Rostoker [16] (KKR), but in a screened representation 
and valid for overlapping MT potentials to leading order in the potential 
overlap. 

Whereas the screened spherical wave must join smoothly onto all passive 
partial waves, we can require that it vanishes at the hard spheres for all the 
active channels except the eigenchannel. This confinement is what makes the 
screened spherical wave localized, provided that localized solutions exist for 
the actual potential, energy, hard spheres, and chosen partition into active 
and passive channels. Since the screened spherical wave is required to vanish 
merely in the other active channels, but not in the eigenchannel, it is an 
impurity solution for the hard-sphere solid and is given by Eq. (JHJ) in TB 
theory. 

Finally, the set of iVMTOs is formed as a superposition of the kinked- 
partial-wave sets for the energies, ei, ejv : 

*52l (r) = S E toBi* r) L { n N l TtnMn ,. (12) 

n=0 Rlfh 

Note that the size of this iVMTO basis set is given by the number of active 
channels and is independent of the number, N + 1, of energy points. The co- 
efficient matrices, in equation (jTJj) are determined by the condition that 
the set of iVMTOs span the solutions, (£j, r) , of Schrodinger's equation 
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with an error 



*W(r) -¥,(£<,!■) = cW(e i -e )(e i -e 1 )...(e i -e N ) (13) 

+o ((e< - e ) (£j - ei) ... (e< - e#)) . 

This condition leads to Eq. (|1U|) and the iVMTO set is a polynomial approx- 
imation for the Hilbert space of Schrodinger solutions, with L„ being the 
coefficients in the corresponding Lagrange interpolation formula. Expres- 
sion (jllj) is for N=l. An iVMTO with N=0 is a kinked partial wave, but an 
iVMTO with N > has no kinks, but merely discontinuities in the (2N + l) st 
radial derivatives at the hard spheres for the active channels. The prefactor, 

, in expression (fT5|) is related to this, [6] and it decreases with the size of 
the set, i.e. with the number of active channels. 

A basis set which contains as many orbitals as there are bands to be de- 
scribed, we shall call truly minimal. For an isolated set of bands, the truly 
minimal iVMTO basis converges to the exact Hilbert space as the energy 
mesh which spans the range of the bands becomes finer and finer. Symmetri- 
cal orthonormalization of the converged 7VMT0 set therefore yields a set of 
atom-centered Wannier functions which are localized by construction. The 
localization depends on the system and on the choice of downfolding. Had 
we, for instance for benzene, chosen instead of the odd sites, sites 1, 2 and 3 
as active, the corresponding iVMTO Wannier functions would have been less 
localized, the Hilbert space spanned by them would have needed a larger N 
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for convergence, and the construction from this set of the maximally local- 
ized Wannier functions would have required a larger cluster. Nevertheless, 
the calculation could have been done. Similarly, in real systems the choice 
of active channels and their hard-sphere radii -their number and main char- 
acters being fixed by the nature of the band to be described- influences the 
properties of the individual iVMTOs, but not the Hilbert space they converge 
to. 

3 Results and Discussion 
3.1 Graphite: a semi- metal 

The bonding in graphite is understood: within a single graphene sheet the 
s, p x and p y orbitals on the two carbon atoms per primitive cell hybridize 
to form a set of sp 2 a-bonding bands which are occupied, and a set of a- 
antibonding bands which are empty. There are two sheets per cell. 

The p z orbitals form a group of bonding and antibonding 7r-bands which 
just touch at the Fermi-level, making graphite a semi-metal. In order to de- 
scribe the set of 7r-bands, we would need to construct two equivalent iVMTOs 
per sheet, each centered on a single carbon atom. We leave it to the method 
to shape the orbitals, subject to the aforementioned boundary conditions 
for the screened spherical waves. The energy mesh must be chosen in such a 
way so that it spans the energy range of the 7r-bands and excludes the energy 
range where the n bands hybridize with other bands. 
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In Figure [T] the band structure of graphite calculated with a full spd 
basis set on each carbon atom is given in black. It is in excellent agreement 
with previous calculations. [17] The red bands have been calculated with 
a basis set comprised of one p z Quadratic-MTO, or QMTO (N—2). The 
energy meshes used for both calculations are given to the right of the band 
structures. The two sets of bands are almost identical, on the scale of the 
figure, with the exception of a small bump at the top of the red bands, where 
hybridization with other bands occurs. Thus, it is possible to describe the 
set of occupied and unoccupied 7r-bands in graphite via just one orbital on 
every carbon atom, shown to the right of the band structure. The orbital is 
localized because it is not allowed to have any p z character on any of the other 
carbons. It is allowed to have other orbital characters, (E.g. s,p x ,d xy ), on 
the other symmetry-equivalent carbons, but such characters are not visible 
in the figure. 

It is even possible to generate orbitals for just the occupied or unoccupied 
7r bands in graphite. In this case we only want to pick half of the bands, and 
therefore we need a basis with, say, a p z orbital on every second carbon atom 
with all other partial waves being downfolded, i.e. passive. Moreover, an 
energy mesh spanning the bonding (anti-bonding) bands must be used in 
order to obtain the bonding (anti-bonding) tc orbital for graphite. Thus, the 
choice of the energy mesh determines which set, bonding or anti-bonding, 
is chosen. In Figures ED and El the orbital on the central carbon atom is 
shown, along with the band structures computed with a full spd basis (in 
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Figure 1: The band structure of graphite calculated with a full spd basis 
is given in black. The red bands have been calculated with a p z QMTO 
(shown on the right) on every carbon atom. The energy meshes used for 
each calculation are given to the right of the band structure. In all figures, 
red denotes a positive and blue a negative isosurface value. For the orbital 
plots, the isosurface values are given in units of a B 3 ^ 2 , where as is the Bohr 
radius. 

black) and those with the truly-minimal basis set we have just specified (in 
red). The agreement between the two sets of bands is excellent, with only 
minor deviations in the upper regions of the downfolded band structure of the 
anti-bonding bands where hybridization with the next higher bands occurs. 
Inspection of the 7r-bonding and anti-bonding orbitals shows that they spread 
out onto the first nearest-neighbour carbon atoms (passive), but they are 
confined not to have any p z character on those carbon atoms, e.g. the second 
nearest neighbours, where the basis set has orbitals (active partial waves). 
The third nearest neighbour atoms also have all partial waves downfolded 
and some p z character may be seen. Clearly, this choice of orbitals breaks 
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Graphite 7t-bonding bands 




Figure 2: As in the previous figure, but the red bands have now been calcu- 
lated with a p z Cubic-MTO (shown on the right) on every second atom within 
a single graphene sheet. The energy mesh is chosen within the occupied part 
of the 7r-band, which is therefore selected. 



Graphite n anti-bonding bands 




Figure 3: As in the previous figure, but the energy mesh is now in the empty 
part of the 7r-band. 

the symmetry; the same Hilbert space would have been obtained had the 
orbitals been placed on the other half of the carbon atoms. 

It is also possible to describe the sp 2 -bonding bands in graphite by placing 
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Carbon s-orbital Carbon p x -orbital Carbon p y -orbital 

g£ gg 85 

isosurface = +/- 0.04 isosurface = +/- 0.04 isosurface = +/- 0.04 



Graphite o-bands 




Figure 4: As in the previous figures, but the red bands have now been cal- 
culated with an s,p x and p y QMTO on every second carbon atom (shown 
above the band structure). Also shown is one of the three congruent sp 2 - 
bond orbitals which arise by symmetrical orthonormalization of the s,p x and 
p y orbitals. 

an s, p x and p y orbital on every second carbon atom, downfolding all other 
channels and using an energy mesh which spans the energy range of the 
bands of interest. The band structure obtained with this basis is given in 
red in Figure 0] and is identical, on the scale of the figure, to the black bands 
which have been calculated using a full spd basis set on every carbon atom. 
The orbitals may spread out onto the nearest neighbour carbons, however 
are confined not to have any s, p x or p y character on the second nearest 
neighbours. Symmetrical orthonormalization of these three orbitals gives the 
well known bond orbital, the carbon sp 2 iVMTO, also shown in the figure. 



19 



The above examples show that it is possible to describe a chosen band, 
or set of bands, with a truly minimal basis set consisting of one iVMTO per 
band. Moreover, the iVMTOs obtained from our method are in-line with an 
intuitive chemical picture of bonding in the solid state, except that they may 
break the symmetry. This arbitrariness originating in the constraint that the 
iVMTOs be atom-centered can be removed by forming linear combinations 
of maximally localized Wannier functions. Hence, iVMTOs should be useful 
not only, for example, as basis sets in linear scaling methods, but also in 
gaining a chemical understanding of periodic systems. 

3.2 Boron Nitride: an insulator 

The bonding in boron nitride is similar to that in graphite: within a single 
layer the s, p x and p y boron and nitrogen orbitals hybridize to form sp 2 
o"-bonding and -antibonding bands. The alternation, however, makes the 
system insulating with a band gap between the bonding and anti-bonding 
7r-bands. In order to describe the occupied bands, it is possible to generate 
either boron or nitrogen centered iVMTOs. It can be expected that the 
electron density, and hence the orbital at a given isosurface, should have a 
maximum closer to the more electronegative element, nitrogen. The method 
needs to do less 'work' if the orbitals are placed initially where the electrons 
are thought to be. Thus, we first place all the orbitals on nitrogen, and let 
the method shape them accordingly. This choice of atom-centered orbitals 
corresponds to the extreme ionic limit, a B 3+ N 3_ configuration. The bonding 
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a and 7r-bands and their respective iVMTOs are shown in the top and bottom 
part of FigureEJ respectively. The s,p x and p y iVMTOs are not shown, since 
they are very similar to those obtained for graphite. The full band structure 
is in excellent agreement with previous calculations. [18] 



BN a- and jr-bands 




Figure 5: The band structure of boron nitride calculated with a full spd basis 
is given in black. The red cx-bonding bands in the upper panel have been 
calculated with an s,p x and p y QMTO on all nitrogen atoms, and the red 
7r-bonding bands in the bottom panel with a p z orbital on all nitrogen atoms. 
Also given is one of the three equivalent sp 2 -bond orbitals and the p z LMTO. 
Boron atoms are purple, nitrogen black. 

In new materials where the bonding is not well understood, it may be 
difficult to decide where the orbitals should be placed. In the following we 
will show that iVMTOs are forgiving: even a bad starting guess can yield the 
correct bands and Hilbert space. Placing all the orbitals on the boron atoms 
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BN a- and rc-bands 




r KM TA L H A 




r KM r A LH A 

Figure 6: As in the previous figure, but with boron-centered a and tc basis 
sets. 

(a B 5 ~N 5+ configuration) yields the bands and orbitals shown in Figure |H1 
The sp 2 -bond orbital looks identical to the one shown in Figured as it should 
be when the energy mesh is converged. For the 7r-bands, more energy points 
are necessary since the orbital has to spread out from a central boron onto 
three neighboring nitrogens. Inspection of the boron (nitrogen) centered 7r- 
iVMTOs makes it plausible that when squared and summed over all boron 
(nitrogen) sites, they give the same electron density, with maxima shifted 
towards nitrogen. 
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3.3 Cesium: a metal 

Within this section we will demonstrate that it is even possible to generate 
Wannier-like functions that span only the occupied bands of a metal. For 
lack of space, we must leave it to the reader to generalize our TB model for 
benzene to the corresponding model for the infinite, slightly dimerized chain, 
and then to study what happens for vanishing dimerization. 

We shall first look at the convergence of the orbital with respect to the 
size of the supercell used. Whereas we can only hope to reproduce the long- 
ranged Friedel oscillations for supercells so large that the facets of the folded- 
in Brillouin zone resemble those of the Fermi surface, much smaller cells turn 
out to reproduce the rough shape of the occupied orbital. This is a manifes- 
tation of what Walter Kohn named the "nearsightedness" of the electronic 
structure of matter. [19] We shall specifically consider cesium because of the 
interesting chemical aspect that its s-electron transforms into a rf-electron 
under hydrostatic pressure. 

The band structure of cesium at ambient conditions (Cs-I) calculated 
with a full sd basis set on every atom is given in Figure [H in black. Superim- 
posed on it in red is the band calculated with one s orbital on every second 
cesium atom, obtained specifically by breaking the symmetry and treating 
Cs metal as CsCl-structured Cs + Cs~. In most regions of the Brillouin zone 
the agreement between the two is good, except it is obvious that with this 
supercell it is not possible to describe the occupied part of the upper band 
along the TX-line. Nonetheless, the orbital can be plotted and is shown in 
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bcc Cs: 2 atoms/unit cell 




Figure 7: Band structure of bcc cesium folded into the 2 atom/cell (CsCl) 
simple cubic Brillouin zone. The black bands were calculated with a basis 
containing the s and d LMTOs on all atoms. The red bands were calculated 
with an s LMTO on every second cesium atom, which is also shown. The 
white atoms have an s orbital placed on them (active), whereas on the black 
atoms all partial waves are downfolded (passive). 

bcc Cs: 16 atoms/unit cell 




Figure 8: As in the previous figure, but for a 16-atom supercell. 
Figure 

The result obtained by doubling the cell in all three directions is shown in 
Figure |H1 Now the occupied band structure has improved considerably, but 
is not yet perfect. The body of the orbital obtained from this calculation is, 
however, very similar to the one generated from the CsCl supercell. The long- 
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ranged tail we do not monitor with the contour chosen in the figures. Already, 
this orbital at low isosurfaces shows the onset of srf-hybridization. At high 
isosurfaces (not shown), the orbital is completely s-like. The hybridization is 
evident in the d z i character found on the nearest-neighbour atoms which have 
all partial waves downfolded, and therefore can posses any orbital character. 
It is a result of the fact that even though the occupied band has primarily s 
character, near the Fermi level some regions with t 2g and d z i character can 
be found. 

Under pressure, cesium undergoes a variety of interesting structural phase 
transitions. At 2.3 GPa Cs-I (bcc cesium) transforms to an fee phase. [20] 
Until recently, it was believed that Cs-II undergoes an isostructural transition 
to Cs-III (which is found in a very narrow pressure range between ~4.2 
and ~4.3 GPa). However, experiments have shown that Cs-III has a very 
complicated structure which is orthorhombic (space group C222i with 84 
atoms per unit cell). [21] At ~4.3GPa Cs-III transforms to the non-close- 
packed tetragonal Cs-IV [22] which at ~12GPa undergoes a transition to 
orthorhombic Cs-V [23] and finally to the double hexagonal close packed Cs- 
VI at ~70GPa. [24] These structural transitions are beleived to be driven 
by the pressure-induced s — > d valence electron transition. [25] Here we will 
generate iVMTOs for Cs-II and Cs-IV in order to visualize how the valence 
s orbital in cesium changes with increased pressure. 

The occupied bands for Cs-II with v/v = 0.6 are still primarily s-like, 
however a fair amount of e 9 -character can also be found. Upon increasing 
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Figure 9: As in the previous figure, but for fee cesium with v/vq = 0.4, and 
a 32-atom supercell. 

pressure, a Lifshitz transition occurs (change in the topology of the Fermi 
surface) as higher lying t 2 g bands cross the Fermi level. Figure |H1 shows 
the band structure of Cs-II for a 32-atom supercell compressed past the 
calculated volume of the Lifshitz transition. The agreement between the 
downfolded and full band structures is excellent and the downfolded bands 
are even able to describe the aforementioned t<i g bands. The orbital at high 
isosurfaces is clearly no longer s-like. It must be noted that the orbital 
calculated for v/vq = 0.6, using the same supercell and downfolding scheme, 
yielded a very similar orbital (not shown), the only difference being that the 
four lobes evident in Figure 01 are located in the same plane. 

Cs-IV with each cesium atom having a coordination of 8 is no longer 
a close-packed structure. It can be viewed as a stacking of prisms with 
a ninety degree rotation from layer to layer in the c-direction. The TB- 
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LMTO calculated charge density fFigurelTUj) shows maxima in the interstitial 
regions, in the center of these prisms. The LMTO obtained by placing one s- 
orbital on every second cesium atom and downfolding all other partial waves 
is given in Figure EH Clearly this orbital can be obtained from that shown 
in Figure |H1 by raising two and lowering two of the lobes. Placing the Cs-IV 
LMTO on all of the active sites and squaring it yields a charge density which 
is almost identical to that calculated with TB-LMTO, [26] giving further 
validation that our method works. Thus we have shown that iVMTOs may 
be used to give a chemical picture of the pressure induced electronic phase 
transition in cesium yielding results which are in-line with those obtained 
from standard electronic structure calculations. 



Cs IV orbital 
isosurface = ♦/- 0.03 



Cs IV electron density 





Cs IV electron density 



Calculated by placing a Cs IV orbital 
on every while atom and squaring it. obtained usinc^TB-LMTO. 




Figure 10: The LMTO obtained for Cs-IV (v/vq = 0.316) with 2-atoms/cell. 
An s orbital was placed on all white atoms (active), whereas the partial waves 
on the black atoms were downfolded (passive). Also shown is the charge 
density obtained by placing the LMTO on every white atom and squaring it, 
along with the charge density obtained from a TB-LMTO calculation. The 
latter charge density is colored by the ELF. The isosurface taken is 0.0063 
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4 Summary and conclusions 



Within this article we have shown that the 3 generation MTO method 
can be used to design a basis set of atom-centered localized orbitals, which 
span the wave functions in a given energy range. For an isolated set of bands, 
arbitrary accuracy may be obtained, with only one basis function per electron 
pair (or single electron, in the case of spin-polarized calculations), simply by 
increasing the density of the energy mesh which spans the energy range of 
the band in question. This method may be applied to insulating and even 
metallic systems to generate Wannier-like functions which are in-line with a 
chemical understanding of bonding in the solid-state. It should therefore be 
a useful analysis tool in, for example, explaining experimental trends for a 
given set of compounds via orbital based arguments; clarifying the bonding 
in novel or even amorphous materials; visualizing pressure-driven electronic 
transitions. 

3 rd generation MTOs may also be useful in generating truly-minimal basis 
sets for order-N methods and in constructing many-electron wave functions 
which can be applied to study strongly correlated systems realistically. In our 
implementation, iVMTOs are generated using the self-consistent potential 
from an LMTO calculation. However, our method may be interfaced with the 
results of any other program, providing that the potential can be expressed in 
terms of a superposition of spherically symmetric potential wells with radial 
overlaps of up to ~60 percent. 
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5 Computational Methods 

Graphite has a hexagonal unit cell. The space group is P6 3 /mmc (194) and 
the two basis carbon atoms are located in the 2b and 2c Wyckoff positions. 
The lattice constants, a and c were taken from experimental data [27] as 
being 2.4642 A and 6.7114 A, respectively. In addition to the atoms on 
the two crystallographic positions, it was necessary to insert two interstitial 
spheres to represent the charge density in the calculation. 

Boron nitride is also found with space group P63/mmc (194), with the 
boron atom located in the 2c and the nitrogen atom in the 2d Wykoff posi- 
tions. The lattice constants, a and c were taken from experimental data [28] 
as being 2.50399 A and 6.6612 A, respectively. It was only necessary to insert 
one interstitial sphere. 

At ambient conditions, cesium crystalizes in the bcc structure. Calcu- 
lations were performed using the experimental lattice parameter of 6.048 
A. [29] Cs-II is fee and the bands shown here were calculated for v/v — 0.4. 
Both Cs-I and Cs-II are close-packed so it was not necessary to insert any 
empty spheres. Cs-IV has the space group I4i/amd with 2 atoms per unit 
cell [22] and it was necessary to insert four empty spheres. 

All of the TB-LMTO calculations [26] were performed using the Vosko- 
Wilk-Nusair (VWN) LDA [30] along with the Perdew-Wang GGA. [31] Scalar 
relativistic effects were included. For graphite and boron nitride a basis set 
consisting of spd LMTOs on the carbon, boron and nitrogen atoms with sp 
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LMTOs on the empty spheres, was used. For cesium, the basis consisted 
of sd LMTOs with pf LMTOs being downfolded. For graphite and boron 
nitride, the calculations utilized 1953 irreducible points in the tetrahedron 
fc-space integrations, [32] 1661 and 897 points were used in the calculation 
for cesium in the bcc and fee structures (one atom per unit cell), respectively, 
and 693 points were used in the calculation for Cs-IV. 

The present version of the iVMTO program is not self-consistent and 
requires the output of the self-consistent potential from the TB-LMTO pro- 
gram. The downfolded band structures are compared with bands computed 
employing a full iVMTO basis set; not with those obtained using the TB- 
LMTO program. In all cases, the default values for the hard-sphere radii, an, 
were used. Thus, all of the hard spheres were slightly smaller than touching. 
However in general, the a R should be taken as 0.9 times the tabulated cova- 
lent, atomic or ionic radii. In the calculations all partial-waves on the empty 
spheres were downfolded. The other downfolding schemes and energy meshes 
employed for particular calculations are given in the results and discussion 
section of the paper. Note that a mesh employing two (N+l=2) energy 
points yields a linear-MTO or LMTO, one with three points a quadratic- 
MTO or QMTO, and so on. More details about the iVMTO formalism can 
be found in [7-9] and references within. 
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